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Abstract 

We present an analytic ansatz to find the effective electrostatic potential and 
Coulomb correlations in multicenter problems, specifically homogeneous and 
doped clusters of metal atoms. The approach is based on a quasi-classical 
density-functional treatment. We focus on the interpretive aspect of our 
findings, particularly on extracting insight regarding the geometric effects 
of Coulomb correlations for any given spatial disposition of ionic cores. For 
singly-doped metallic clusters we obtain a direct visualization of the varia- 
tions of both screening and Coulomb correlations with changes of location of 
the dopant atom. This analysis provides a way to interpret recent observa- 
tions of the variability of physical properties of metal clusters with changes 
of composition and geometry. 
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Collective effects induced by Coulomb correlations in atoms have been studied in two 
ways. In the first, both hydrodynamic theory and local approximate dielectric theory have 
been used; neither of these takes into account either shell structure or the single-particle 
spectrum of the valence electrons. JTJ These methods are capable, at most, of giving gross 
trends in dynamical properties. The second route instead uses a fully quantal description 
based on the one-electron excitation spectrum and corresponding wavef unctions. A recent 
collection of papers provides a description of methods and results of the application of 
many-body techniques in atomic theory. || 

The way electrons are correlated can be inferred from the probability distribution implied 
by their wavefunction. To make it possible to make such inferences, however, we must be 
a bit thoughtful about how we present this distribution. Even for a two-electron atom, we 
begin with a function of six independent variables in a fixed center-of-mass system. We would 
like to extract from this a description in no more than two or three independent variables, 
something we can represent pictorially and visualize. For a three-body system such as He** 
or the valence electrons of Mg, a natural and practical way to carry out such a description 
has emerged as an analytic reduction of the probability density | \I/ (r 1; r 2 ) | 2 to the joint 
probability density p (n, ?"2, O12), where 612 is the angle between the position vectors T\ and 
r 2 of two electrons. This in turn makes it straightforward to compute and display the 
conditional probability density d (r 2 , 612', ri), for finding one of the two electrons at distance 
r\ from the nucleus, and at an angle 6 12 from the vector from the nucleus to electron 2 if 
electron 2 is at a distance r 2 from the nucleus. Three-dimensional graphs of d (r 1; r 2 , #12) 
provide a vivid and precise way to depict the correlation of two electrons. H With this 
probability density d, one can compare wave functions of different qualities, see what roles 
long-range and short range correlation play in various states, exhibit the relative importance 
of angular and radial correlation, and compare correlation in different atomic systems. The 
work, started in the late 1970s by Rehmus, Kellman, Roothaan and Berry |3],|J, provides a 
generalization of other quantitative descriptions of electron correlation. f^-Q 

For the system of more particles, we have yet to find a comparably powerful approach 
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because so much information is contained in the wave function and we do not know how to 
extract what is relevant in a manner adaptable to pictures. An exception emerges, obviously, 
for the high-density limit of the Fermi fluid where a collective description of electrons is 
likely to be optimal. This collective description is based on the organized behavior of the 
electrons brought about by their long-range Coulomb interactions. The long-range Coulomb 
interactions, subject to the screening among the electric charges, act to couple together the 



motion of many electrons, giving rise to the well-known quantum density oscillations. |L5 
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Working toward extracting insights from the probability density and making use 
of the quasi-classical description for the (valence) electron gas, we develop here an analyt- 
ical ansatz which allows us to find and visualize the effective electrostatic potential and 
Coulomb correlations in multicenter problems. We apply this ansatz to the case of mod- 
erately large metallic clusters. To anticipate what follows, let us state our findings: by 
using a generalized partition function for valence electrons (the Bloch density matrix), the 
electron self-distribution in the common potential V (r) is derived in terms of many-body 
perturbation theory. |T(J This approach produces the electron density p (r) as a functional 
of V (r) (with p(r) ~ V(r)), which is valid for describing metallic systems, i.e. systems 
with a high-density valence electron gas. Further, inside the electron gas of density p(r), 
we introduce the cluster cage formed by the positive ion cores with the spatial distribution 
given by p + (Hi), and apply Poisson's equation to the cluster as a whole. (The vectors Rj 
are position vectors of the ions.) The self-consistent solution of this equation gives the col- 
lective description of the cluster constituents, electrons plus ions. This generalization of the 
Coulomb interaction results in a superposition of quantum oscillations given by long-range 
contributions and screening on the smooth " semi classical" potential ||17|| . We focus on their 
interpretive aspects and specifically on extracting insights regarding the geometric effects 
of Coulomb correlations for any given spatial disposition of ionic cores. Also, we explore 
the case of a foreign metal atom doping a otherwise-homogeneous cluster of metal atoms. 
The approach presented here provides us with a direct visualization of the way both the 
screening effect and the Coulomb correlations change with changes of the location of the 



impurity. This analysis is important in the context of recent observations of the role played 
by composition and geometry in changing the physical properties of metallic clusters fl8|-|2"0|| . 

Consider a fixed positive ion distribution in space p + (Hi), with Rj the positions of the 
ions measured from the center of the cluster. The ion cluster cage has a net charge measured 
in appropriate units equal to zN, where z is the electric charge of one ion (for simplicity, we 
shall restrict the discussion to single-valent metals, z — 1) and N denotes the total number 
of ions in the cluster. If a gas of valence electrons carrying an equal number of negative 
charges is introduced, so that the system is strictly neutral, then the electrons redistribute 
themselves so as to shield the positive charges at large distances and minimize the Coulomb 
self-energy of that gas, and also satisfy the Fermion constraints on the electrons. In the 
high density limit of the Fermi fluid, Bohr's Correspondence Principle applies and we can 
introduce the collective description of the electron gas based on the methods of statistical 
mechanics. The generalized partition function of the valence electrons moving in the common 
potential V (r) can be written in terms of the wave functions (r) and energy levels £, as 

r(r',r,/3) = ^**(r')^(r) e -^ , (1) 

where (3 = (ksT) 1 , ks is Boltzmann's constant and T the absolute temperature. By 
integrating this along the diagonal, where r' = r, we obtain the ordinary partition function 
of statistical mechanics. Eq. (1) is the Bloch density matrix and if we operate with the 
one-particle Hamiltonian 

H s = - 1 -V 2 r + V(r) , (2) 

on T, and compare the result with that obtained by differentiating T with respect to (3, then 
we find the Bloch equation 

".r = -g, (3) 

which has the form of the time-dependent Schrodinger eq., with f3 playing the role of it. 
The boundary condition required to define the solution of (3) follows from the completeness 
theorem for eigenf unctions, namely 



r (r', r,0) = 6 (r' - r) . 

In the high density limit, the behavior of the electrons is simple, and the Coulomb interaction 
can be treated as a perturbation of the motion of the free electrons. Therefore the solution 
of eq. (3) may be written as 

r(r',r,/3) =r (r',r,/3)- (4) 
J dr" f Q dp'T (r', r",/?-/?') V (r") T (r", r,/3') , 

where T is the Bloch density matrix for an assembly of free electrons. The derivation of 



Dirac's density matrix 7 (r', r) from (4) is described in Ref. [16]. It consists of an iterative 
procedure to find the perturbation terms in T. The first step is replacement of To for F in 
the integral and integration over f3' . If the Bloch density matrix V is determined, Dirac's 
density matrix 7 (which is actually the electron density) may be obtained by using the 
Laplace transform relation connecting V (r', r,(3) and 7 (r', r,£) 



/■OO 

r(r',r,/3) =(3 / #7 (r', r,£) exp 



where £ is an intensive energy variable conjugate to (3. 

Following the above procedure, the electron density is obtained as 

p{r)=p Q -—J dv V{v) — - 2 . (5) 

fc 3 

Here kp is the Fermi wavevector, p is the free-particle density, po = 3^, and ji (x) is the 
first-order spherical Bessel function. 

With the primary form (5) , the electron density p (r) makes further mathematical com- 
putations very difficult. To go further, we need to simplify by linearizing. Usually, this 
linearization proceeds by adding the assumption that V (r) varies slowly in space, the 
Thomas- Fermi approximation. Accordingly, V (r') is replaced in (5) by V (r) [^TJ. One 
obvious point needs to be stressed here: the discrete positive ion distribution used here 
produces a Coulomb potential with a far more rapid spatial variation than that of the 



5 



frequently-invoked continuum distribution of the "jellium" models p2| . The linearization in 
(5) can still be made under the following assumptions. Our main observation at the outset is 
that in a high-density electron gas, any electric charge is screened out very rapidly, namely, 
at distances beyond a characteristic Debye screening length, say q^ 1 (which is inversely pro- 
portional to \/Fp, as can be seen below). Also, we notice that the slow variation of V (r) is 
usually supposed to be over a de Broglie wavelength for an electron at the Fermi surface, 
that is 27r /kp. From this view, we may say that a possible conflict with the use of the 
Thomas-Fermi approximation occurs only close to the positive ions, closer than a shielding 
distance . Disregarding this limitation, we use the Thomas- Fermi approximation in (5) 
to obtain, after a straightforward integration, 

P (t)= Po -^V{t) , (6) 

where q% = ^£., with a# the Bohr radius. The results are not strongly affected by this 
approximation. For example, good agreement within natural limits has been obtained pre- 
viously |2Bj for the fullerene molecule described in this way and without the simplifying 



linearization. At the same time, we can see that the density follows the potential closely, 
which means that the validity of the theory is ensured, as we already stipulated above, by an 
appropriate requirement on the electron density. Of course this theoretical model loses its 
validity at large distances from the ion locations because the electron density vanishes, and 
at very short distances, towards the center of the cluster cage, where the density becomes 
infinite with the potential |24[ . 



Eq. (6) is a quasi-classical result obtained in the high density approximation for the 
valence electrons. Within the quasi-classical approximation, |25] local variations of the 



electron density leave the exchange contribution unchanged, as a consequence of its nonlocal 
(quantum) character. Therefore, a first-order quantum correction to this quasi-classical 
result represents the exchange energy 
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Consequently, we may assume that the common potential V (r) is generated only by the 
electron distribution in the presence of the discrete ionic background. We may therefore set 
up the basic Poisson equation to yield 

N 

Ay = 47rp -g V(r)-47r5>5(r-R i ) . (7) 

i 

We have to solve a self-consistent field problem which accounts for the electron distribution 
profile in the presence of a discrete positive background. The last term on the right side of 
eq.(7) represents the density of positive charge with Ri the average distance of an ion from 
the center of the cluster and % is an index running over the ions, each with electric charge 
z%. These locations are chosen without regard to the stability of the configuration. 
According to the principle of superposition, Poisson's equation (7) may separate: 

N 

A{V 1 + V 2 )=±irp -q 2 V 1 -q 2 V 2 -±ir'£z i 6{r-R i ) , (8) 

i=i 

which means that we have to solve two simpler equations rather than one very complex 
equation. The first is given by 

AV 1 = Anp - q 2 V x , (9) 

and represents the effective electrostatic potential due to the electron self-distribution where 
the discrete nature of the positive charges is disregarded. This equation will be solved inside 



a large sphere of radius R, which has to contain most of the valence electron density. f26 
The second equation becomes 

N 

AV 2 = -q 2 V 2 -A7iJ2z l 5(r-R i ) , (10) 

i=i 

and accounts for the remaining terms of the total potential. The discrete nature of the 
positive background is employed here. 

By Fourier transformation, the latter equation becomes 

r 1 N r 

J dk V 2 (k) (k 2 - q 2 ) exp (zkr) = — £ z { J dk exp [ik (r - R*)] , (11) 

i = l 



wherefrom 



1 " exp(-ikRi 



2**U 1 W-ql 

and the potential is simply 

1 * r exp[ik(r-R,)] 



or 

V 2 (r) =-J2ziJ2i l I dk-^j^klr-Rtl) 
tr J k 2 -q 2 

J dn k Y t * m (9 k , ip k ) Y lm (9i, ifi) , 

in terms of spherical Bessel functions ji (k |r — Rj|). After the integration over fl k , the above 
equation for the potential reduces to 

2 N r k 2 

V 2 (r) = -5>/dA;- ^sin^lr-R,!) . (13) 

vr ^ J k 2 -q 2 

The former equation (9) (subject to appropriate boundary conditions, as we will see 
below) deals with the effective Coulomb potential due to the electron distribution in the 
super-sphere of effective radius R. To solve it we exploit the fact that (9) separates in 
spherical polar coordinates r, 6, if. The solution of Poisson's equation (9) is given by 

V 1 (r,9, ( p) = '£F lm (r)Y lm (e, ( p) . (14) 

l,m 

Each Fi m is actually independent of m and satisfies the radial equation 

;£<riS)-^iS=4*A.-iF. ■ (15) 

Hence we now drop the subscript m. Strictly, terms corresponding to / = have been 
considered separately in solving the above equation. The general solution for the radial 
equation is 

rp ( \ (47r) 3/2 A 00 . / N , -Boo / \ / ic x 

Fi{r) = 5 — PoH sin (g r)H cos (q r) (16) 

q r r 

1 C ■ 

+ J2 i * J , i B i™ ex P (-^or) + (-I) 1 A lm exp (q r) 

j=0 q r3+ 



where 

/(/ + l)(Z+j)! 

C "' _ 2^! (/-,)! ' (17) 
and A o, -Boo; ^Zm and B im are constants that will be determined. The effective Coulomb 
potential V\ can be written then as 

T r 47T A 00 B 00 

Vi = -yPo H sin (q r) H cos (q r) (18) 

qz r r 

1 C ■ 

+ E ' E ~T~fel exp (-g r) + A /m exp (g r)] F im (0, y?) , 

where the prime in the right hand term of the equation means that the summation over / 
begins from / = 1. This potential has to be finite for r = 0, which means that B 00 = 0, and 
Bim = (— l) l+1 Ai m . 

Taking (12) into account, we find the total effective Coulomb potential inside the super- 
sphere is 

V in = % Po + — sin (q r) + V 2 (r) (19) 
% r 

+ E ' E Clj j { ~^ A im [exp (q r) + exp (-q r)] Y lm (9, cp) , 

lm j=0 % r 

everywhere except for r = Rj. Outside the super-sphere, a Laplace equation applies and the 
solution vanishing at infinity is 

K* = ^ + E'^(M • (20) 

' lm 1 

If the potential is specified on the surface of the bounding sphere, the coefficients entering 
(17) and (18) can be determined by evaluating V (R, 9, ip) and using 

A lm = J dnY lm (9 : <p)g(9 : <p) , (21) 

where g (9, ip) is an arbitrary function. Here, g represents a "pseudo-charge density" designed 
to be a smooth, nodeless function which, in order to maintain the electrical neutrality of the 
entire system, has to agree exactly with the true charge density outside the region bounded 
by the super-sphere of radius R. 
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An additional comment is appropriate here regarding the present theory. We begin by 
asking, "How unique is the potential in eqs. (18) and (19)?" If we demand that our "pseudo- 
charge density" g agrees with the true charge density outside the "super-sphere" then the 
potential is uniquely determined. The inside region is not uniquely fixed by this procedure; 
however, if we require that the total charge of the valence electrons be normalized, then the 
fraction of electronic charge contained in this region must be large, e.g., more than ~ 95% 
of the total. This means that the behavior of the charge contained in this region must 
dominate the static properties of the metallic cluster. 

The effective cluster potential given by (18) displays the usual collective aspects of the 
electron gas. The primary manifestations of the collective behavior are a) collective oscil- 
lations of the valence electrons as an entity, the so-called "plasma" oscillations, and b) the 
screening of the field of any individual electric charge beyond a characteristic length g^ 1 . 
The former is fundamentally a diffraction effect, the electron wave nature being essentially 
disregarded in this kind of calculation. The screening of the ionic fields causes the remainder 
of the electron gas to stay diffuse, and so leads to a deficiency of negative charge just out- 
side the immediate neighborhood of each positive ion enclosed in its neutralizing, co-moving 
electronic cloud. Thus, the cluster potential exhibits additional spatial oscillations which 
are not determined solely by its behavior in the neighborhood of r. In a collective oscil- 
lation, each individual electron suffers a small periodic perturbation of its velocity (recall 
that the electron density (5) is a result of the perturbation of the kinetic operator) and 
position due to the combined potential of all the other particles, both positive and negative. 
The cumulative potential of all the electrons may be quite large since the long range of the 
Coulomb interaction permits all the electrons to contribute to the potential at every point. 
The collective behavior of the electron gas dominates phenomena involving distances greater 
than the characteristic length q^ 1 , while the individual particle component is associated with 
the random thermal motion of the electrons. In the approximate level of this analysis, the 
effects of collective excitation on the correlation are neglected, as a second-order effect. 

Usually, the long range of the Coulomb interactions having the character of (18) precludes 
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immediate application of these results to the calculation of the ground-state energy of the 
cluster. Therefore, we are not able to perform a minimization of the ground-state energy 
with respect to the volume of the super-sphere. Consequently, the self-consistency of the 
potential is affected by this lack of information. 

With all the assumptions of the model and its mathematical output now presented, 
we may already point out some general characteristics we may expect for the behavior of 
the effective cluster potential. Since in the present perturbation approach the expansion 
of the electron density is based on plane waves, the cluster potential displays a high value 
in the central region. The potential is strongly dependent on the Coulomb correlations 
and, naturally enough, very sensitive to the position of the positive ions. (We discuss this 
aspect later.) Nonlocal effects due to the particle spins in the mean field for electrons are 
disregarded. Hence the method produces state-independent potentials. 

In the following discussion, we work out an example of a metallic cluster M± 3 with 
icosahedral symmetry that closely approximates spherical symmetry. The model for the 
ionic cores is that of hard-spheres occupying a total volume in space equal to Vt ions . The 
valence electrons are highly confined between the ionic cores. The unperturbed density p 
is expressed by 

3 N 

*±/i / s 3 L ions 

where r s is a point in the space available to the electrons, the "electronic interspace" , outside 
the ion cores. This means that we have subtracted from the entire volume of the super- 
sphere of radius R the volume assigned to the ionic cores VL ions ] N is the total number of 
the delocalized electrons, equal to the number of ionic charges. The distance between the 
centers of the central and outer ions is the bond length. For numerical calculation we set 
Ri = 5.9 a.u.. The core volume VL ions is usually computed by taking into account the ionic 
radius; In our first example, we let this radius be 2.7Aa.u. so the resulting volume VL ions 
is 1122(a.-u.) 3 . For the other two cases, this volume is a free parameter. The electronic 
interspace was chosen to be r s = 0.75 a.u., in accordance with the high-density electron 
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gas requirement (r s <C 1), and by imposing that 95% of the total electrons must be inside 
the super-sphere, the super-sphere radius becomes R = 6.5 a.u. In Fig. 1 we can see the 
corresponding effective potential inside the super-sphere as a function of r and 9. In Fig. 
2, the spatial dependence of the same potential is displayed along the coordinates 9 and 
at the radius where the outer ions lie, Ri = 5.9 a.u.. These pictures show the regular, 
collective characteristics we discussed above. The oscillations we observe are a manifestation 
of the self-consistency of Poisson's equation and represents the main correlation effect of the 
electron gas in the metallic state (the high-density limit). 

Despite the collective aspects which contribute to the mean-field character of the effective 
potential, V (r) remains sensitive to the geometry and composition of the cluster. If, for 
example, one host atom in the cluster cage is replaced by an impurity atom A, the potential 
reflects this structural change. We explore this property in the following. Let us assume that 
the impurity A is a trivalent metal atom (A) which releases all three of its valence electrons 
into the Fermi sea in the bulk volume. We take its core volume to be the same as the M 
ions. The substitution does not change the symmetry and we assume that the most stable 
geometry of the cluster has the dopant at the centre. Figs. 3 and 4 show the corresponding 
effective electrostatic potential for the AM U system. By comparing this with the effective 
mean-field potential for M13, we observe that the presence of the trivalent atom at the 
center makes the potential much deeper. (See Fig. 3, where V is displayed as a function of 
r and 9 .). The delocalized electrons polarize inward, toward the high Coulomb field of the 
central, trivalent ion. Fig. 4 shows also that the shape of the quantum oscillations at the 
cluster surface changes relative to the previous case. The amplitudes of oscillations become 
rather uniform which means that the screening among the electric charges is much better 
for this system than for the homogeneous cluster, a consequence of the increased number 
of delocalized electrons, from 13 to 15. Moreover the better screening effect here results in 
a change in the effective force acting in the electron gas. This can be seen in Fig. 4 as a 
phase-shift of electron density oscillations at the position of the positive ions. (Compare 
this with the effective potential displayed in Fig. 2.) 
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If we move the trivalent ion to the outer shell, the disturbance of course goes toward 
the surface (see Fig. 5) and the electron density is enhanced around the vertex where 
the trivalent ion is located. This behavior of the effective potential is supplemented by 
the appearance of more pronounced Coulombic correlations of the valence electrons near the 
surface. The quantum oscillations are sensibly disturbed by the trivalent impurity located on 
the cluster surface. This disturbance appears as irregular behavior along the 9 coordinate at 
constant r. Also, a large potential difference, about 1 a.u., can be seen in Fig. 6 between the 
position of the trivalent impurity (8 = 0, ip = 0) and the antipodal position (6 = tt, ip = 0) 
occupied by a host ion. The potential difference leads to a displacement of the electronic 
cloud towards the position of the trivalent ion and a deficiency of negative charge in the 
opposite direction. Consequently, a diffusive trend of electron density oscillations can be 
observed in the hemisphere at 9 = n. 

We may conclude that the effective electrostatic potentials for metallic clusters is subject 
to important Coulombic correlation effects which can be visualized at the proper scale by 
employing a discrete description for the positive background. The electron density shows 
a static screening which is rather localized near the positive charges and supplemented by 
the long-range oscillatory behavior. For singly-doped binary metallic systems, the depth of 
the effective electrostatic potential depends on the dopant position in the cluster geometry. 
The collective aspects of the excitations of electrons delocalized through the cluster volume 
are strongly perturbed by the presence of the impurity. The main effect of these structural 
rearrangements of the ions is the change of the effective potential, as we have shown. The 
change of the cluster potential will alter, in turn, the ordering of the related electron shell, 
a fact which has been observed in many experiments [|18| -|20||. Therefore our findings may 
be interpreted as qualitative support for various models explaining the shell inversions for 
doped metallic clusters [p7]-|3"0|j. Obviously other kinds of changes of dopant atoms may 
induce still different effects, that will depend on their locations in the cluster. 

Finally, we may say that the method developed here is simple and flexible and can 
yield, to some extent, accurate approximations to the exact effective potentials with minor 
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computing effort. Also, it has the advantage of physical immediacy, i.e., the present approach 
is easy to interpret. This makes the method useful for a fast check of the effective potential 
to systems, clusters of heavy elements, for example, presently beyond the capability of more 
accurate approaches. 
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Figure captions 



Fig. 1 

The effective electrostatic potential inside the cluster cage for M i3 for < r < 6.5 a.u., 
< 9 < 7r rad and <p — 0. 
Fig. 2 

The same effective electrostatic potential displayed as a 9 — p> plot at the position of the 
surface ions. 
Fig. 3 

The corresponding effective electrostatic potential for AM 12 system for < r < 6.5 a.u., 
< 9 < Ti rad and ip — with the impurity in the center of the icosahedral cluster cage. 
Fig. 4 

The effective electrostatic potential showed in Fig. 3 displayed here along the coordinates 
9 and ip at the position of surface ions. 
Fig. 5 

A plot analogous to that of Fig. 3, for the AM\ 2 system with impurity at the vertex, 
[9 = 0, ip = 0) and for < r < 6.5 a.u.. 
Fig. 6 

The 9 — <p spatial dependence of the potential displayed in Fig. 5 at the position of the 
outermost ion shell. The electrons must be in their ground state and highly confined in the 
cluster cage. 
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